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We introduce two sandpile models which show the same 
behavior of real sandpiles, that is, an almost self-organized 
critical behavior for small systems and the dominance of large 
avalanches as the system size increases. The systems be- 
come fully self-organized critical as the system parameters are 
changed, showing that these systems can make a bridge be- 
tween the well known theoretical and numerical results and 
what is observed in real experiments. A simple mechanism 
determines the boundary where self-organized criticality can 
or cannot be found, which is the existence of local chaos. 



The concept of Self-Organized Criticality (SOC) was 
introduced by Bak, Tang and Wiesenfeld (BTW) in 1987 
to denote a phenomenon in which out of equilibrium, 
multidimensional systems, drive themselves to a critical 
state characterized by a power-law distribution of event 
sizes [jl| . Until then, the studies of fractal structures were 
related to equilibrium systems where this kind of distri- 
bution appears only at special parameter values where a 
phase transition takes place. In that pioneering work, the 
concept of SOC was illustrated by a model for sandpiles 
and since then an enormous amount of research on SOC, 
both theoretically and experimentally, has been done. 
Among other phenomena in which SOC has been con- 
nected with are earthquakes M and evolution S. 

The existence of SOC in an experiment with a quasi- 
one-dimensional pile of rice was demonstrated by Frette 
et. ed 0]. They found that the occurrence of SOC de- 
pends on the shape of the rice. Only with sufficient 
elongated grains avalanches with a power-law distribu- 
tion occurred. For more symmetric grains a stretched 
exponential distribution was seen. Christensen et. al || 
introduced a model for the elongated rice pile experiment 
in which the local critical slope varies randomly between 
1 and 2. They found that their model, known as the 
Oslo rice pile model, reproduced well the experimental 
results on the quasi-one-dimcnsional rice pile. In a recent 
publication ||, we introduced a fully deterministic one- 
dimensional SOC system, which presents the same qual- 
itative and quantitative behavior of the Oslo system. In 
other words, they belong to the same universality class. 

When one goes to sandpiles with geometry of two- 
dimensions a different picture emerges. That is, the 
models pi predict the presence of power-law distribu- 
tions and the experiments do not display them. The 
most well known sandpile experiments can be classified 
in two types: (a) local dropping of sand in the center 



of the pile (jj and (b) uniform driving, more specifi- 
cally the rotating drum experiments [@|§. In class (a) 
it was found that small systems present scaling almost 
consistent with SOC, but in large systems another regime 
with big avalanches belonging to a different distribution 
appears. In class (b) one sees small avalanches that 
seems to display a power-law distribution of limited size. 
These small avalanches are interwoven with big system- 
wide avalanches which belong to a different distribution. 
There are studies of avalanches in natural settings pOl , 
which also show similar distribution as the ones observed 
in the sandpile experiments. 

Although SOC has not been observed in sandpile ex- 
periments is well known that power law distributions does 
exist in Nature, one of the most well known cases being 
Gutenberg-Richter law for earthquakes [0| in agreement 
with SOC. It is clear that there is a missing piece in 
this puzzle. That is, would there be simple models that 
would display the observations in real sandpiles, and still 
present SOC in other parameter regions? To our knowl- 
edge it still missing in the literature such models and it 
is the aim of this letter to introduce them. Our models 
makes a bridge between the SOC and non SOC behavior, 
and the boundary that characterizes the separation be- 
tween the two regimes is the absence or presence of local 
chaos. However, we will see that local chaos is a necessary 
but not sufficient condition for SOC to be destroyed. 

Our models are inspired by the model introduced in 
Ref. H , and is also governed by coupled lattice map (that 
is, systems characterized by discrete time and continu- 
ous values for the space variables). Here we increase the 
dimensionality and change the drive and the relaxation 
rules. Having in mind the sandpiles experiments we first 
introduce a model for the local dropping of sand. We 
assume a two-dimensional square lattice of linear size L 
and to each site i,j in the lattice there is associated to 
it a variable Xij with x G [0, +oo), which is to represent 
the local slope of the pile. The dynamics of the model is 
described by the following algorithm: 

(1) Start the system by assigning random initial values 
for the variables Xjj , so the they are below a chosen, 
fixed, threshold x t h- 

(2) Choose a nearly central site of the lattice and up- 
date it slope according to x%j — Xth- 

(3) Check the slope in each element. If an element i,j 
has Xij > Xth, update x%j according to x! i A = 4>(xij — 
Xth), where <f> is a given nonlinear function that has two 
parameters a and d. Increase the slope in all its nearest 
neighboring element according to x' nn = x nn + Ax/4, 



where Ax 



x\ ■ and nn denotes nearest neighbors. 



(4) If x'ij < x t h for all the elements, go to step (2) (the 
event, or avalanche, has finished). Otherwise, go to step 
(3) (the event is still evolving). 

Without losing generality, we can take Xth = 1- In 
our simulation in step (2) we have chosen the site with 
i = j = L/2. The nonlinear function we use is 



4>{x) 



1 — d — ax, if x < (1 — d)/a, 
0, otherwise. 



(1) 



The parameter d would be associated with the minimum 
drop in energy after an event involving one single el- 
ement and a would be the parameter associated with 
the amount of friction between the grains. That is, the 
smaller the a, the larger the friction and the smaller the 
change in the slope of the pile. We have tested several 
other functions and found that the quantitative and qual- 
itative results we show here are robust. The important 
ingredient being the shape of tfi(x) in the vicinity of x = 0. 
In contrast with the one-dimensional case Q, it is not 
required here that (j)(x) be periodic in order to find the 
presence SOC. 

We have chosen to evolve the system using parallel dy- 
namics with open boundary conditions. It is beyond the 
scope of the present letter to study the several possible 
variations of our models. Further results on these models 
will be presented in a future publication Il2]. The dis- 
tribution of time duration of the avalanches will also be 
presented in the future, but our preliminary results show 
that they are qualitatively similar to the ones for the size 
distribution. 

We display an example of our simulations in Fig. |lj, 
where we show the distribution of events P(s) involving 
s update steps, that is, the size of the avalanche. The 
events that involve all the elements of the system have 
been excluded from our analysis. However, in this model 
of local dropping of sand we have observed that noth- 
ing very distinct will occur if they are also included in 
the statistics. In (a) we show P(s) for L = 64, d = 0.1 
and vary a and in (b), we use a = 1.5 and vary d, keep- 
ing L — 64. We notice the existence of two regimes. 
For small d, power-law distributions, that is SOC, ap- 
pears only when a < 1, whereas if d is large, we see SOC 
even with a > 1. As in real sandpiles M, the almost 
SOC regime is characterized by a region with an appar- 
ent scaling for small events, and the big events belong 
to a different distribution. To illustrate this, we show 
in Fig. RYc) simulations with a = 1.5 and d — 0.1 and 
varying L. In small systems P(s) can be fit to a scaling 
form of the type 



P{s 1 L)=L- l 'G{s/L v ), 



(2) 



as shown in Fig. H(d), where we have used /? = 6 and 
v = 3. The function G is not well fit by a power-law, 
since one can clearly see in the figure that it is curved. 
We have found fO that it is consistent with a stretched 



exponential, as in real sandpiles |l3[ . The observations 
of Himalayan avalanches [|l0| also display the kind of dis- 
tribution shown in Fig. UKc) for large systems. 

The second model we introduce here is for the rotating 
drum experiment [pip), which we call uniform drive, since 
the slope of the pile increases uniformly for all the grains. 
The algorithm is similar to the one described above, with 
the exception of step (2), which is now replaced by 

(2) Find the element in the lattice that has the largest 
x, denoted here by x max . Then update all the lattice 
elements according to Xij = x-ij + x t h — x max . 

We show examples of P(s) for this model in Fig. 0. 
There, in (a) we fix d and vary a, and in (b) we fix a and 
vary d. In both cases we have used L = 64 and the events 
that involve all the elements of the system have been 
excluded. Distinctly from the model of local dropping, it 
seems here that there is a power-law distribution for any 
parameter value. However, the behavior is not exactly 
SOC. We have found that SOC is only seen if a < 1 or 
d is large enough, as in the case of the local dropping. 
When a > 1 and d is smaller than a given value, we see a 
SOC-like behavior only for small values of L. As L grows, 
there is a transition to a different behavior, in which the 
larger the system, the smaller the power-law region, as 
Fig. g(c) shows. 

System-wide avalanches have been reported in the ro- 
tating drum experiment [pLp| that belong to a different 
distribution than the one of the small avalanches. This is 
exactly what we see in this model for a > 1 and small d. 
In Fig. 0(a) we show all the events of the system including 
the ones that involve all the elements (for a — 3, d = 0.1 
and varying L). We see that the scaling region does not 
get bigger as the system size increases, and a peak related 
to the events involving all the elements is seen. We have 
found that the intermediate size events (10 < s < 100) 
can be fit by a scaling of the type P(s, L) — L~ l3 G{s/ L v ) 
with (3—1 and v = 0, as shown in Fig. J3|(b). As in the 
rotating drum experiment |9|,|l3| the function G in this 
case is closer to a power-law than in the case of the ex- 
periment with local dropping of sand. 

For given a and d the slope of the power-law distribu- 
tion seems to be the same for both models, that is, local 
dropping of sand and uniform drive, but the slope varies 
with a, as shown in Fig. fil and 0. Consequently, the uni- 
versality class of these models vary with the parameters. 
In Fig. ||(a) we show the scaling given by Eq. (||) with 
(3 = 3.55 and v — 2.70 for the local dropping (Id) and 
/3 = 3.55 v = 2.85 for the uniform drive (ud). In the limit 
a — ► or d — * 1 the slope of the power-law distribution 
is the same as the one in the BTW model m and the 
conservative OFC model 0, that is, P(s) ~ s^" 1245 . 

We next investigate what would happen if the relax- 
ation function is just a random number generator. In 
other words, instead of using (j) in step 3 of the above al- 
gorithm we now use x[ = p, where p is a random number 
uniformly distributed in the interval [0, 1 — d]. We have 



found SOC for any d e (0, 1] with the same exponents 
as the BTW model. This is displayed in Fig. [|(b) where 
we show the size distribution for d = 0.01 in the cases 
of local dropping (Id) and uniform drive (Id) . Therefore 
nonlincarities in <f> and consequently non ergodicity are 
necessary for the SOC behavior to be destroyed in these 
models. In that figure, we also show the case in which 
x\ = 0, which corresponds to d = 1. In this limit we 
recover the OFC M model for the case of uniform drive. 

The reason why a = 1 determines a special boundary, 
in which SOC may or may not be present, is due to the 
fact that it marks the boundary in which local chaos ex- 
ists. By "local" here we mean that is is at the grain level, 
no matter what happens at the system level. Chaos is 
defined as the exponential divergence of trajectories that 
start with almost the same initial conditions. In fact, 
if we consider two copies of a given system, copy 1 and 
copy 2, with all the elements having the same x except 
that in copy 1 the element x\a = X*, whereas in copy 2 
x? = x* + A. It is not difficult to find that after one 
iteration by <j), the separation of the two elements instead 
of A will be aA. So, if a > 1 the separation increases 
(that is, we have local chaos) whereas if a < 1 the separa- 
tion decreases. We have seen from the above results that, 
since SOC can exist even with a > 1, if d is large enough 
or L is smaller than a given value, this means that local 
chaos is a necessary but not sufficient condition for SOC 
to be broken. Consequently, the true boundary between 
SOC and non SOC depends on three parameters, that is, 
a, d and L. Further studies are necessary to understand 
this boundary more clearly. 

In summary, we have introduced a model for sand- 
piles, or other systems that present avalanche like be- 
havior, that reproduce very well the observation in real 
sandpiles. We have found that a simple mechanism, i.e 
local chaos, can explain the breakdown of SOC in those 
systems. Based on our results, we believe that with an 
appropriate choice of grains SOC will be seen even in real 
sandpiles. Grains with large friction would be the best 
candidates for this. 
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FIG. 1. (a) Frequency of events involving s updates for 
(a) variable a with d = 0.1, and (b) variable d and a = 1.5 
with L — 64. In (c) we vary the system size and use a = 1.5 
and d — 0.1, and in (d) we show the fitting using the scaling of 
Eq. 2 for small systems, with (3 = 6 and v = 3. In this figure 
we are using the model for local dropping of sand. We have 
used 4 x 10 6 avalanches in all the simulations of this paper. 



FIG. 2. (a) Frequency of events involving s updates for 
(a) variable a with d = 0.1, and (b) variable d and a = 1.5 
with L — 64. In (c) we vary the system size and use a = 1.5 
and d — 0.1. In this figure we are using the model for uniform 
drive (that is, the rotating drum experiment). 



FIG. 3. (a) Frequency of events involving s updates for 
a — 3 and d — 0.1 with varying L. The peaks in the distribu- 
tion are the events involving all the elements of the system. 
In (b) we show the fitting using the scaling of Eq. 2 for the 
region of avalanches of intermediate sizes, using j3 — 1 and 
v = 0. In this figure we are using the model for uniform 
drive. 



FIG. 4. (a) Frequency of events involving s updates for 
a = 0.3 and d — 0.1 with varying L using the scaling of Eq. 
2. For the case of local dropping we have used /? = 3.55 
and v = 2.75. For uniform drive we have used j3 — 3.55 and 
v = 2.85. In (b) we show the cases in which after relaxation 
the variable x is a given random number uniformly distributed 
between and 0.99 for local dropping of sand (Id) and for 
uniform drive (ud). We also show the cases in which x is 
relaxed to after an event, which is the same as having d = 1. 
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